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Part  I 


Computational  Modeling  of  the  VKI 
Inductively  Coupled  Plasma  (ICP)  Facility 
with  the  Minnesota  US3D  Code 

1  Motivation 

he  US3D  code  developed  at  the  University  of  Minnesota  in  Prof.  Graham  Candlers  research 
group  has  gained  widespread  use  throughout  NASA  and  U.S.  Air  Force  research  laboratories. 
It  is  a  state-of-the-art  flow  solver  for  challenging  hypersonic  flows  and  has  been  validated 
with  high  enthalpy  shock  tunnel  experiments  (CUBRC  and  Caltech  tunnels  for  example) 
as  well  as  with  flight  data  for  over  a  decade.  This  work  has  been  extremely  successful  and 
the  next  major  challenge  is  to  develop  new  finite  rate  gas-surface  interaction  models  and 
incorporate  them  into  the  US3D  code.  This  would  enable  high-fidelity  solutions  to  flows 
involving  surface  catalysis  and  ablation;  problems  of  increasing  importance  to  the  USAF. 
Recently,  Prof.  Thomas  Schwartzent rubers  group  at  the  University  of  Minnesota  implemented 
a  Finite- Rate-Catalytic  (FRC)  boundary  condition  into  the  US3D  code.  This  is  a  high-fidelity 
gas-surface  interaction  model  and  the  goal  is  to  determine  the  model  parameters  and  validate 
the  US3D  solutions  with  high  quality  data  from  VKI  experiments. 


Gas-surface  interactions  are  difficult  to  study  in  shock  tunnels.  Although  the  flow  is  character¬ 
istic  of  flight  conditions,  the  test  time  is  on  the  order  of  microseconds  and  material  surfaces  do 
not  have  time  to  heat  up  as  they  would  during  flight.  In  the  VKI  Inductively  Coupled  Plasma 
(ICP)  facility,  the  flow  is  subsonic  and  as  a  result  the  test  times  are  such  that  the  material 
reaches  temperatures  expected  in  flight,  and  high-temperature  gas-surface  interactions  (cataly¬ 
sis  and  ablation)  can  be  carefully  observed  over  long  test  times.  If  the  boundary  layer  conditions 
in  the  subsonic  ICP  flow  can  be  matched  to  the  boundary  layer  conditions  in  hypersonic  flow  at 
flight  conditions,  then  the  expected  behavior  of  the  TPS  during  flight  can  be  inferred  directly 
from  ICP  experiments.  Specifically,  VKI  employs  the  Local  Heat  Transfer  Simulation  (LHTS) 
concept  to  match  boundary  layer  edge  properties  in  an  ICP  experiment  to  the  corresponding 
hypersonic  flight  conditions. 

There  are  two  main  goals  to  this  Minnesota- VKI  collaborative  project: 

•  To  use  the  US3D  code  to  validate  the  LHTS  approach  employed  by  VKI.  Since  LHTS  is  a 
very  useful  and  widely  used  technique  used  to  interpret  ICP  data,  a  better  understanding 
of  its  accuracy  and  potential  limitations  would  be  a  valuable  contribution. 

•  To  use  the  US3D  code  to  directly  interpret  ICP  data  and  thereby  develop  new  high- 
fidelity  models  for  high-temperature  gas-surface  interactions.  Basically,  in  the  same  way 
that  US3D  has  been  rigorously  validated  with  shock-tunnel  data  (for  developing  gas- 
phase  models),  we  plan  to  do  the  same  with  ICP  data  for  gas-surface  interaction  models. 
If  possible,  such  general  models  could  then  be  applied  using  US3D  to  a  large  number  of 
hypersonic  problems. 
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As  a  reference,  an  example  flow  field  solution  obtained  from  the  US3D  code  (nonequilibrium), 
overlaid  with  a  local  thermodynamic  equilibrium  (LTE)  solution  from  VKI  for  the  case  of  a 
plasma  jet  flowing  over  a  probe-geometry  is  shown  below  in  Fig.l. 


Figure  1:  Representativecomparison  between  nonequilibrium  (US3D)  and  equilibrium  (VKI- 
LTE)  solutions  for  an  ICP  flow  field. 


2  Project  Accomplishments 

The  ICP  flow  is  subsonic  and  thus  the  plasma  jet  is  expected  to  be  very  close  to  thermal 
and  chemical  equilibrium.  The  important  non-equilibrium  physics  thus  only  occur  within  the 
boundary  layer.  Thus  in  all  simulations,  the  focus  is  on  the  boundary  layer  profiles  and  the 
boundary  layer  edge  conditions.  Indeed  the  LHTS  technique  is  used  to  determine  the  boundary 
layer  edge  conditions  from  an  ICP  experiment,  after  which  hypersonic  flight  conditions  are 
determined  that  result  in  precisely  the  same  boundary  layer  edge  properties.  Specifically,  the 
boundary  layer  edge  conditions  along  the  stagnation  streamline  are  defined  by  several  key 
parameters;  namely,  total  enthalpy  (H),  temperature  (T),  bulk  flow  velocity  (u),  and  the 
velocity  gradient,  ,  where  is  the  radial  velocity.  In  fact,  the  LHTS  defines  the  boundary  layer 
edge  as  the  point  of  inflection  in  the  profile  (seen  in  upcoming  figures).  Thus,  the  simulation 
results  presented  in  this  section  generally  focus  on  profiles  of  these  four  variables  along  the 
stagnation  line  approaching  the  test  article. 

The  US3D  code  has  not  been  validated  or  tested  for  subsonic  flow  simulations  prior  to 
this  research.  The  numerical  methods  used  by  US3D  are  in  fact  optimized  specifically  for 
hypersonic  flow  and  thus  much  of  the  research  discussed  in  this  report  is  preliminary  work 
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aimed  at  methodically  validating  US3D  for  subsonic  plasma  flows. 


2.1  Grid  convergence  and  US3D  numerics  for  simple  subsonic  flows 

A  grid  convergence  study  was  performed.  A  typical  grid  density  commonly  employed  for  VKI 
LTE  solutions  is  shown  in  Fig. 2  below.  US3D  simulations  were  performed  using  this  grid  and 
two  refined  grids  at  twice  and  four  times  the  resolution.  The  resulting  profiles  of  u,  T,  H, 
and  are  plotted  in  Fig. 2.  There  is  a  noticeable  difference  in  stagnation  line  profiles  between 
the  baseline  grid  and  the  two  refined  grids.  As  a  result,  it  is  recommend  to  increase  the  grid 
resolution  by  roughly  a  factor  of  two  for  future  ICP  simulations. 


Figure  2:  Grid  convergence  study  on  stagnation  line  profiles  (left). 

Baseline  mesh  density  (right). 

Furthermore,  a  series  of  US3D  simulations  were  performed  to  assess  general  convergence  be¬ 
havior,  influence  of  domain  size,  and  influence  of  the  numerical  stencil  order-of-accuracy,  on 
simple  subsonic  flow  solutions.  As  seen  in  Fig. 3,  a  uniform  high-temperature  air  flow  (typi¬ 
cal  of  ICP  conditions:  T=5840  K,  u=209  m/s,  p=2000  Pa,  and  Twall=300  K)  was  simulated 
over  10cm  diameter  sphere.  It  was  determined  that  given  uniform  far-held  subsonic  conditions, 
that  the  US3D  code  was  able  to  rapidly  converged  to  a  solution  with  no  modifications  at  all. 
The  solution  was  independent  of  the  domain  size  (domains  of  10,  20,  and  40  sphere  diameters 
were  tested).  Finally,  it  was  found  that  various  US3D  stencils  resulted  in  essentially  the  same 
solution.  For  example,  a  6th-order  accurate  stencil  did  not  noticeably  change  the  solution  com¬ 
pared  to  more  standard  2nd-order  accurate  stencils.  This  is  because  the  higher-order  stencils 
available  in  US3D  are  designed  for  supersonic  turbulent  hows  where  steep  gradients  perme¬ 
ate  the  how.  However,  this  subsonic  plasma  jet  how  is  very  smooth  in  comparison  and  so  we 
conclude  that  the  use  of  a  2nd  order  stencil  is  accurate.  Finally,  it  is  interesting  to  note  that 
for  the  simulations  depicted  in  Fig. 3,  there  was  unsteady  behavior  in  the  wake  of  the  sphere. 
Although  this  was  simulated  only  approximately  by  US3D,  the  wake  had  no  noticeable  effect 
on  the  stagnation  line  prohles  (the  forebody  how). 
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Figure  3:  US3D  numerics  study  for  simple  subsonic  flow  over  a  sphere. 


Conclusions: 

•  The  baseline  VKI  mesh  is  too  coarse,  we  recommend  using  2x  the  resolution 
(  4x  the  number  of  cells  for  2D/Axi  simulations). 

•  The  unsteady  wake  behind  a  sphere  (only  approximately  modeled  by  US3D), 
was  shown  to  have  no  influence  on  the  stagnation  line  profile.  This  supports 
the  use  of  steady-state  solutions  for  these  flow  problems. 

•  It  was  determined  that  the  grid  resolution  in  the  boundary  layer  is  very  im¬ 
portant  for  solution  accuracy,  however,  using  6th  order  interpolation  schemes 
provides  no  additional  benefit  over  standard  2nd  order  methods  for  this  type 
of  flow. 

2.2  Comparison  of  2D  vs.  3D  simulations 

The  VKI-LTE  simulations  are  axis-symmetric  simulations.  However,  some  experimental  sample 
geometries  of  interest  are  not  axis-symmetric.  For  example,  a  flat  plate  geometry  with  a  blunted 
leading  edge  has  been  tested  in  the  VKI  ICP  facility  to  study  off-stagnation  heating  and  surface 
catalysis.  This  flow  is  actually  three-dimensional,  however,  VKI  currently  assumes  that  a  2D 
flow  solution  is  accurate.  Since  the  US3D  code  is  a  highly-scalable  3D  code,  we  ran  both  2D 
and  3D  simulations  with  US3D  and  compared  the  boundary  layer  results.  The  3D  simulation 
required  over  3  million  elements  and  the  solution  was  obtained  using  256  core  CPUs  (solution 
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seen  in  Fig. 4  below).  The  stagnation  line  profiles  were  found  to  be  virtually  identical  between 
2D  and  3D  solutions,  as  was  the  solution  on  the  plane  of  symmetry. 


I 

$ 

U 

t 


Figure  4:  Stagnation  line  profiles  (u,  T  and  dv/dx)  compared  for  2D  and  3D  simulations. 


Conclusions: 

•  The  boundary  layer  profiles  along  the  stagnation  line  and  the  symmetry-plane 
solutions  are  virtually  identical  between  2D  and  3D  simulations,  supporting 
the  use  of  the  2D  simulations  used  by  VKI. 

•  However,  these  US3D  simulations  did  not  include  a  round  plasma  jet  imping¬ 
ing  on  the  3D  plate  geometry.  This  would  add  another  aspect  to  the  3D  flow 
and  these  simulations  should  be  repeated  with  a  more  realistic  plasma  jet  to 
be  sure  the  2D  simulations  are  still  accurate. 

2.3  Subsonic  Outflow  and  Plasma  Jet  Inflow  Boundary  Conditions 

To  accurately  model  the  ICP  flow  field,  a  realistic  plasma  jet  profile  must  be  specified  as 
the  inflow  boundary  conditions  to  a  US3D  simulation.  Likewise,  proper  subsonic  boundary 
conditions  must  be  specified  within  US3D  in  order  to  accurately  model  the  ICP  flow. 


A  new  subsonic  outflow  boundary  condition  was  implemented  and  tested  in  the  US3D  code. 
The  boundary  condition  extrapolates  two  characteristic  variables  from  the  interior  solution 
to  each  outflow  cell  face.  These  two  characteristics  together  with  a  prescribed  background 
pressure  are  then  sufficient  to  determine  the  complete  gas  state  and  these  state  values  are 
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imposed  in  ghost-cells  bordering  the  outflow  cell  faces. 


Next,  the  US3D  code  was  modified  to  specify  an  arbitrary  plasma  jet  profile  as  an  inflow 
boundary  condition.  The  LTE  code  developed  at  VKI  is  capable  of  modeling  the  joule-heating 
inside  the  inductively  coupled  plasma  torch;  basically  the  process  that  generates  the  plasma. 
This  code  has  been  developed  over  a  number  of  years  and  validated  with  data  from  the  VKI  ICP 
facility.  It  is  not  necessary  or  desirable  for  US3D  to  have  this  joule-heating  modeling  capability. 
Since  the  resulting  plasma  exits  through  a  jet  at  equilibrium  conditions,  this  is  a  more  relevant 
inflow  condition  to  prescribe  as  a  US3D  boundary  condition.  At  some  position  in  the  jet  flow 
(outside  of  the  torch),  the  plasma  jet  profiles  are  extracted  from  the  VKI  LTE  solution.  It  is 
noted  that  in  the  future  it  may  be  possible  to  experimentally  measure  the  plasma  jet  properties 
and  directly  use  these  as  inflow  conditions.  This  simulation  setup  was  portrayed  earlier  in  Fig.l. 


Subroutines  were  written  to  extract  the  plasma  properties  at  any  x-location  in  the  jet  from 
a  VKI  LTE  solution  and  automatically  prescribe  these  as  inflow  conditions  within  a  US3D 
simulation.  Since  the  VKI  LTE  code  and  US3D  code  do  not  necessary  employ  the  same 
equilibrium  constants  in  their  chemistry  models,  tests  were  performed  to  ensure  that  the 
equilibrium  profiles  predicted  by  LTE  were  the  same  as  predicted  by  US3D. 


An  example  solution  corresponding  to  an  11  species  air  model  is  shown  below  in  Fig.  5.  Here 
the  solid  lines  are  the  mass  fractions  extracted  from  the  VKI  LTE  solution  at  a  specific  x- 
location.  Therefore,  these  solid  lines  represent  the  inflow  boundary  conditions  for  the  US3D 
simulation.  The  symbols  show  the  final  converged  US3D  solution  for  all  mass  fractions  at  the 
same  x-location.  Since  the  mass  fractions  remain  virtually  unchanged,  this  demonstrates  that 
the  initial  conditions  extracted  from  the  VKI  LTE  solution  are  consistent  with  the  equilibrium 
predictions  of  the  US3D  code.  The  progression  from  pure  O  and  N  in  the  jet  core  (small  y)  to 
pure  02  and  N2  in  the  ambient  air  far  from  the  jet  (large  y)  is  evident  from  Fig. 5. 

Conclusions: 

•  A  new  characteristic  subsonic  boundary  condition  was  implemented  in  the 
US3D  code. 

•  A  new  plasma  jet  profile  inflow  boundary  condition  was  implemented  in  the 
US3D  code.  Equilibrium  chemistry  profiles  were  verified  to  be  in  agreement 
with  VKI  LTE  solutions. 

2.4  VKI  LTE  vs  US3D  solution  comparisons 

The  main  aspect  of  comparing  VKI  LTE  solutions  to  US3D  solutions  is  that  US3D  is 
a  nonequilibrium  code  that  does  not  employ  a  local  thermodynamic  equilibrium  (LTE) 
assumption  and  also  simulates  the  flow  through  the  boundary  layer  edge  to  the  surface  of 
the  object,  and  therefore  does  not  couple  to  a  boundary  layer  code  at  the  boundary  layer 
edge.  Such  a  comparison  could  aid  in  the  validation  of  the  LHTS  technique  employed  by  VKI 
researchers  and  also  ultimately  validate  a  single  consistent  model  (US3D)  for  the  entire  flow 
field  that  can  be  used  to  interpret  the  experimental  data. 
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Figure  5:  Mass  fraction  for  11-species  air.  Lines  are  extracted  from  VKI  LTE  solution  (US3D 
inflow  conditions).  Symbols  are  final  converged  US3D  solution. 


General  contour  comparisons  are  shown  below  in  Fig. 6  for  temperature  and  bulk  velocity.  The 
VKI  LTE  solution  (blue  lines)  extends  through  the  entire  domain  (torch,  jet,  and  test  object), 
whereas  the  US3D  solution  is  initialized  in  the  jet.  In  general,  good  agreement  is  found  for 
contours  of  T  and  u  in  the  jet  and  it  is  expected  to  see  differences  near  the  body  since  one 
solution  is  an  LTE  solution  and  the  US3D  solution  is  nonequilibrium. 


One  interesting  test  case  that  was  investigated  was  performing  US3D  simulations  of  the  same 
flow  with  different  inflow  boundary  locations.  As  seen  in  Fig. 7,  plasma  conditions  were  ex¬ 
tracted  from  the  VKI  LTE  code  at  two  x-locations,  and  used  as  boundary  conditions  for  two 
US3D  simulations.  The  agreement  between  the  two  US3D  results  strongly  supports  the  assump¬ 
tion  that  the  plasma  jet  flow  is  indeed  in  thermal  and  chemical  equilibrium.  Specifically,  the 
contours  between  Cutl  and  Cut2  in  Fig. 7  are  from  a  nonequilibrium  US3D  solution,  yet  as  the 
flow  reaches  the  location  of  Cut2,  the  solution  still  agrees  with  the  LTE  boundary  conditions 
extracted  along  Cut2  from  the  VKI  LTE  code.  Furthermore,  both  US3D  solutions  continue 
to  agree  downstream.  This  result  demonstrates  that  according  to  US3D,  the  flow  is  in  equi¬ 
librium,  otherwise  the  two  US3D  solutions  would  have  started  to  deviate  and  this  deviation 
would  propagate  downstream.  Finally,  the  stagnation  line  profiles  are  shown  in  Fig. 8  for  both 
US3D  solutions  and  they  are  found  to  agree  almost  exactly. 

One  possible  source  of  discrepancy  comes  from  the  top  boundary  condition.  In  reality,  the 
plasma  jet  exits  into  a  large  chamber  of  stagnant  air  which  would  drive  a  large  recirculating 
flow  region  to  develop  between  the  jet  and  the  chamber  walls.  To  include  these  regions  in 
an  unsteady  CFD  solution  would  be  very  difficult  in  terms  of  CPU  time  and  unsteady  CFD 
modeling.  Thus  the  top  boundary  condition  is  set  somewhat  arbitrarily  and  a  small  study 
was  completed  to  investigate  this  issue.  Specifically,  converged  solutions  were  obtained  when  a 
bulk  flow  velocity  of  u=10,  20,  and  80m/s  was  imposed  along  the  top  boundary.  Although  it  is 
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Figure  6:  (top) Temperature  field  contours  (blue-VKI  LTE,  red-US3D). 

(bottom)  Velocity  field  contours  (blue-VKI  LTE,  red-US3D). 


unclear  precisely  what  the  flow  properties  would  be  along  this  top  boundary,  it  is  clear  that  this 
would  be  a  region  of  low  momentum,  low  temperature  flow,  compared  to  the  high  momentum, 
high  temperature  jet  flow,  and  thus  may  not  influence  the  solution.  Indeed,  we  found  that 
the  velocity  imposed  along  the  top  boundary  has  little  to  no  influence  on  the  stagnation  line 
flow  properties  (or  flow  near  the  test  object)  as  long  as  the  imposed  velocity  is  low  (uj80m/s). 
The  higher  the  velocity  set  at  the  top  boundary,  the  faster  the  convergence  of  the  solution,  so 
typically  a  velocity  of  u=20m/s  was  used  for  the  remaining  simulations.  Figure  9(left)  shows 
the  3  jet  profiles  resulting  from  these  three  velocity  boundary  conditions  1cm  upstream  of  the 
test-object.  Figure  9 (right),  shows  the  stagnation  line  profiles  of  the  variables  of  interest,  where 
no  noticeable  difference  is  observed. 

The  most  important  comparison  between  US3D  and  VKI  LTE  solutions  is  again  the  bound¬ 
ary  layer  variables  along  the  stagnation  streamline.  These  are  plotted  in  Fig.  10.  Here  we  see 
excellent  agreement  for  bulk  velocity  and  density,  however,  we  see  large  discrepancies  for  tem¬ 
perature  and  total  enthalpy.  US3D  is  a  nonequilibrium  code,  that  employs  a  two-temperature 
model  (T  and  T^).  We  expected  to  see  thermal  equilibrium  for  this  flow  and  therefore  T=T^ 
at  virtually  all  points  in  the  flow.  However,  we  see  anomalous  behavior  of  T^  along  the  plasma 
jet  core,  near  the  inflow.  Furthermore,  we  see  a  7/discrepancy  in  total  enthalpy  persist  through- 
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Cut  1 


Cut  2 


Figure  7:  Two  US3D  solutions  (velocity  contours)  initialized  at  two  different  locations  in  the 
plasma  jet. 


dv/dr 


Figure  8:  Stagnation  line  solutions  for  the  two  US3D  solutions  shown  in  Fig. 7 


out  the  flow.  This  enthalpy  discrepancy  is  even  evident  in  the  first  few  cells  near  the  inflow 
boundary  condition,  despite  both  simulations  specifying  the  same  plasma  inflow  conditions. 

After  substantial  investigation,  code-to-code  comparisons,  and  modification  of  physical  models, 
it  was  determined  that  all  of  these  discrepancies  are  due  to  the  fact  that  US3D  does  not  cur¬ 
rently  model  electronic  energy.  If  the  US3D  solution  is  post-processed  to  account  for  electronic 
energy  modes,  the  computed  enthalpy  at  the  inflow  moves  into  near  perfect  agreement  with 
the  VKI  LTE  code  as  seen  in  Fig.  11.  However,  this  post-processing  assumes  LTE  at  each  point 
in  the  US3D  solution  to  determine  the  enthalpy.  Since  we  are  attempting  to  prove  the  LTE 
assumption  is  valid,  it  is  required  to  add  a  nonequilibrium  model  for  electronic  energy  to  the 
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Figure  9:  (left)  jet  profile  solution  1  cm  upstream  of  the  test-object.  Three  profiles  are  shown 
corresponding  to  three  top-boundary  velocity  values  (u=  10,  20,  80  m/s). 

(right)  Stagnation  line  profiles  corresponding  to  the  three  solutions. 
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Figure  10:  Stagnation  line  solutions  from  the  VKI  LTE  code  (lines) 
and  from  the  US3D  code  (symbols). 


US3D  code. 

Adding  a  model  for  electronic  energy  would  also  eliminate  the  anomalous  behavior  of  T^.  The 
problem  is  that  in  the  plasma  jet  core  the  temperatures  are  so  high  that  the  air  is  completely 
dissociated  (mass  fractions  of  diatomics  are  close  to  machine  zero).  This  situation  has  never 
been  encountered  by  the  US3D  code  since  in  hypersonic  flows  the  gas  is  never  completely  dis¬ 
sociated  due  to  nonequilibrium  effects.  Numerically,  the  vibrational  temperature  is  computed 
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Figure  11:  Stagnation  line  enthalpy  and  temperature  profiles  for  US3D  with  electronic  energy 
post-processing  compared  to  the  VKI  LTE  solution. 


from  the  vibrational  energy  considering  the  diatomic  species  that  this  energy  is  distributed 
among.  However,  in  the  jet  core  there  is  essentially  no  vibrational  energy  and  more  trouble¬ 
some,  there  are  no  diatomics  present.  As  a  result,  the  value  of  T^  is  not  well  defined  and 
numerically  its  value  has  large  error.  To  some  degree  the  value  of  T^  doesnt  influence  the 
solution  since  there  is  almost  no  vibrational  energy  present  in  that  region  anyway,  however,  the 
problem  is  that  the  dissociation/recombination  model  employed  by  US3D  is  Parks  T-Tv  model. 
So  Tvib  is  specifically  used  to  set  the  dissociation  and  recombination  rates.  Thus  even  if  the 
anomalous  T^  value  does  not  influence  the  vibrational  energy  equation,  it  strongly  influences 
the  chemistry.  As  a  result,  it  is  important  to  modify  US3D  to  avoid  such  numerical  problems. 
By  adding  electronic  energy  and  coupling  this  to  the  vibrational  temperature  this  would  mean 
that  even  with  no  diatomics  there  will  be  electronically  excited  atoms  present.  This  electronic 
energy  would  then  result  in  a  well-defined  value  for  the  combined  vibrational-electronic  tem¬ 
perature. 

Conclusions: 

•  Good  agreement  between  US3D  and  VKI  LTE  solutions  was  found  for  hydro- 
dynamic  properties  of  the  plasma  jet  and  boundary  layer.  However,  significant 
discrepancies  were  found  in  the  translational  and  vibrational  temperatures 
and  total  enthalpy  throughout  the  simulation  domain,  including  near  the  in¬ 
flow  boundary. 

•  It  was  determined  that  these  discrepancies  are  due  to  the  fact  that  US3D 
does  not  currently  model  electronic  energy. 

•  It  was  determined  that  implementing  electronic  energy  and  coupling  with 
vibrational  energy  (a  combined  electronic-vibration  temperature)  will  bring 
the  enthalpies  into  agreement  and  also  eliminate  numerical  issues  regarding 
the  vibrational  temperature. 
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3  Future  Work 


This  collaborative  project  that  involved  a  small  amount  of  initial  funding  to  have  Prof. 
Schwartzentruber  visit  and  work  at  VKI  for  3  months  is  the  start  of  a  much  longer  term  and 
productive  collaboration  between  the  University  of  Minnesota  (computational  expertise)  and 
VKI  (experimental  expertise).  This  research  is  directly  related  to  a  large  Multi-University- 
Research  -Initiative  (MURI)  funded  by  AFOSR  and  based  at  the  University  of  Minnesota  on 
Fundamental  Processes  in  Hypersonic  Flows.  Collaborative  work  will  continue  on  this  specific 
project  over  the  next  year. 


Specifically,  we  are  in  the  process  of  adding  a  model  for  electronic  energy  to  the  US3D  code.  We 
are  also  implementing/testing  a  variety  of  gas-surface  interaction  boundary  conditions.  Within 
the  next  few  months  we  expect  to  have  full  US3D  solutions  of  the  ICP  flows  with  all  required 
physical  models  for  both  gas-phase  and  gas-surface  reactions.  At  this  point,  direct  comparison 
with  ICP  experimental  results  will  be  performed.  We  expect  two  publications  towards  the  end 
of  2012.  The  first  involving  validation  of  the  LHTS  technique  by  full  nonequilibrium  US3D 
simulation  and  the  second  involving  interpretation  of  experimental  measurements  (stagnation 
and  off-stagnation  data)  for  the  creation  of  new  gas-surface  interaction  models. 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 

15 


Part  II 


Numerical  simulation  of  plasma  jet  and 
gas-surface  interaction  around  blunted 
plate  with  COOLFluiD 

Introduction 

This  document  presents  the  computational  model  (physico-chemical  model  and  numerical 
method)  and  some  preliminary  results  for  the  numerical  analysis  of  a  Thermal  Protection 
System  (TPS)  test  in  the  VKI  Plasmatron  wind  tunnel.  The  experiment  consists  of  a  high- 
enthalpy  plasma  jet  impinging  on  a  flat  plate  with  blunted  leading  edge,  featuring  different 
Thermal  Protection  System  materials  with  different  catalytic  properties.  The  numerical  simu¬ 
lation  has  been  performed  with  the  aerothermodynamic  code  inside  the  COOLFluiD  platform. 
Both  nonequilibrium  and  gas-surface  interaction  effects  have  been  taken  into  account  in  our 
investigation. 

0.1  COOLFluiD  platform 

COOLFluiD  (Computational  Object  Oriented  Libraries  for  Fluid  Dynamics)  [1,  2,  3,  4]  is 
VKI  collaborative  software  environment  for  high-performance  scientific  computing  where  dif¬ 
ferent  numerical  techniques,  physical  models,  post-processing  algorithms  can  coexist  and  work 
together.  Herein,  each  numerical  method  or  physical  model  is  encapsulated  into  an  independent 
dynamic  module  (or  plug-in  library)  that  can  be  loaded  on  demand  by  user-defined  applications. 
Some  of  the  main  features  of  COOLFluiD  include: 

•  different  multidimensional  parallel  solvers  for  compressible  and  incompressible  flows, 
based  on  different  space  discretization  techniques  for  unstructured  meshes,  including 
cell-centered  Finite  Volume  (FV),  Residual  Distribution  Schemes  (RDS),  Spectral  Finite 
Volume/Difference  and  standard  Finite  Element  (FE). 

•  explicit  (Runge  Kutta  n-order)  and  implicit  (3  point-backward,  Crank-Nicholson,  time- 
limited)  time  marching  schemes; 

•  interfaces  to  several  linear  system  solver  packages  such  as  PETSc,  Trilinos,  Pardiso,  etc. 

•  compressible  Reynolds  averaged  Navier- Stokes  with  k— cj,  BSL,  SST  and  Spalart-Almaras 
turbulence  models; 

•  multi-temperature  and  Collisional  Radiative  models  for  arbitrary  gas  mixtures  in  thermo¬ 
chemical  nonequilibrium; 

•  Magneto-Hydrodynamics  models  for  Space  Weather  applications; 

•  coupling  algorithms  for  multi-physics  and  multi-domain  simulations,  allowing  both  strong 
and  loose  coupling 

•  robust  mesh  deformation  algorithms. 
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1  Physical  modeling 


The  sets  of  governing  PDE’s  (Navier-Stokes)  describing  a  2D  axisymmetric  flow  in  a  continuum 
regime  can  be  expressed  in  conservative  and  hypervectorial  form  as: 


dU  drF  drF %  drFcr  drFt  drF<t 


+ 


+ 


+ 


+  S 


(1) 


dP  dt  dx  dr  dx  dr 
where  x  and  r  are  the  axial  and  radial  directions.  Herein,  according  to  the  terminology  adopted 
in  COOLFluiD,  U  are  the  conservative  variables,  P  the  update  variables,  Fc  and  Fd  respectively 
the  convective  and  diffusive  fluxes,  S  the  source  term. 


1.1  Park’s  2-temperature  model 

In  the  case  of  a  gas  mixture  in  thermal  and  chemical  nonequilibrium  (TCNEQ),  the  conserva¬ 
tive  and  update  (or  natural)  variables  corresponding  to  a  2-temperature  model  [5,  6]  can  be 
expressed  as  follows: 


u  =  (ps,pu,pE,pev)T,  P  =  (ps,u,T,Tv)T  (2) 

where,  in  particular,  ps  represents  the  partial  densities,  ey  is  the  vibrational-electronic  en¬ 
ergy  per  unit  mass  and  Ty  stands  for  the  vibrational-electronic  temperature.  In  this  case,  the 
convective  and  diffusive  fluxes  are  defined  by: 
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where  the  pressure  p  is  given  by  the  following  equation  of  state: 
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where,  in  particular,  pe  represents  the  electron  pressure,  ys  are  the  species  mass  fractions,  Ms 
the  species  molar  masses  and  Te  =  Ty  if  thermal  nonequilibrium  is  assumed. 

The  tensor  of  viscous  stresses  f  appearing  in  Eq.  3  is  defined  as 
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and  it  is  computed  using  Stokes’  hypothesis  of  negligible  bulk  viscosity  effects.  Herein,  the  dy¬ 
namic  viscosity  fi  is  computed  rigorously  from  kinetic  theory  by  using  the  transport  algorithms 
described  in  [7,  8]. 

The  mass  diffusion  fluxes  ps uf  are  computed  by  solving  the  Stefan-Maxwell  system  [9,  10] 
of  equations  which  consist  of  a  linear  system  (in  the  diffusion  fluxes)  of  as  many  equations  as 
the  chemical  species  are  present  in  the  mixture.  This  system  is  supplemented  by  the  auxiliary 
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condition  that  the  sum  of  the  diffusion  fluxes  is  zero.  Moreover,  by  imposing  the  ambipolar 
constraint  stating  that 


ne  =  rii=  Y  ns  (6) 

s=ions 

where  ne  is  the  electron  number  density  and  rii  is  the  ionic  number  density,  the  flow  field  in  an 
ionized  gas  mixture  can  be  considered  electrically  neutral  [9,  11]. 

The  roto-translational  and  vibrational-electronic  heat  fluxes  q  and  q y  are  defined  as 


AVT  —  Ay  VTy  —  'y  ^  psushs 

s 

(7) 

\vVTv  -^Ps^sK 

(8) 

s 


where  the  A  is  the  roto-translational  thermal  conductivity,  Ay  is  the  vibrational-electronic 
thermal  conductivity  and  the  species  enthalpies  hs  are  given  by 


ha  =  hst{T)  +  hse(Tv)  +  h) 

atoms 

(9) 

hs  =  hst(T)  +  hsr(T)  +  hsv(Tv)  +  hse(Tv)  +  hsf 

molecules 

(10) 

hs  =  hst{Tv) 

free  electrons 

(11) 

where  /if ,/if  ,/if ,  /if,  hSf  are  respectively  the  translational,  rotational,  vibrational,  electronic  and 
formation  enthalpies  for  species  s. 

1.1.1  Mass  production  /  destruction  terms 

The  source  terms  in  Eq.  2  can  be  expressed  as 
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(12) 


In  the  first  contribution,  due  the  axisymmetric  formulation,  an  additional  viscous  stress  com¬ 
ponent  tqq  in  the  circumferential  direction  6  reads: 


2  f  du  dv  v 
^0  = --/W  —  +  —  -  - 


(13) 


The  mass  production/destruction  term  us  for  chemical  species  with  partial  densities  ps  which 
appears  in  Eq.  12  is  formulated  as  follows: 
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where  the  forward  reaction  rates  kfr  corresponding  to  a  5-species  neutral  air  mixture  used  in 
this  work  are  taken  from  [12]. 
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1.1.2  Energy  exchange  terms 


Qvt,  the  energy  exchange  (relaxation)  between  vibrational  and  translational  modes  due  to 
collisions  can  be  expressed  as 


nvt 


(15) 


where  e^C  is  the  equilibrium  vibrational  energies  of  molecules  m  evaluated  at  the  roto- 
translational  temperature,  according  to  the  Landau- Teller  formulation  [13].  The  latter  assumes 
mono-quantum  energy  transfers:  in  the  collision  between  two  molecules  or  between  a  molecule 
and  an  electron,  one  colliding  particle  can  gain  or  lose  only  one  energetic  level,  while  the  energy 
of  the  other  particle  remains  unchanged.  The  relaxation  time  rm  is  given  by  Millikan  and  White 
[14]  with  Park’s  correction  for  high  temperatures  [15]: 


TPark 

Tm  =  r^w(p,T)  +  {amcmnm)~l  (16) 

where  <jm  is  the  effective  cross  section  for  vibrational  relaxation  processes,  cm  is  the  average 
molecular  velocity  of  molecule  m  and  nm  is  the  number  density. 

Finally,  £lcv  stands  for  the  vibrational  energy  lost  or  gained  due  to  molecular  dissociation  or 
recombination: 

V,CV  =  J2“mDm  (17) 

m 

Among  the  several  possibilities  for  the  choice  of  Dm  reported  in  literature,  the  simplest  one  is 
to  impose  Dm  =  e^,  with  being  the  vibrational  energy  of  the  molecule.  Qet  represents  the 
energy  exchange  due  to  inelastic  collisions  between  the  electrons  and  the  heavy  particles: 
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where  z/e,s  is  the  effective  collision  frequency  of  electron  with  heavy  particles  as  defined  in  [5]. 
fl7  corresponds  to  the  energy  loss  due  to  electron  impact  ionization  and  is  given  by 


^  =  (19) 

ions 

where  he?s  is  the  molar  rate  of  production  of  species  s  and  Is  is  the  energy  lost  per  unit  mole 
by  a  free  electron  in  producing  species  s  by  electron  impact  ionization  [11]. 

For  a  broader  introduction  to  the  thermodynamic,  chemical  and  transport  properties  modeling 
employed  in  this  work,  the  reader  may  refer  to  [5,  4,  7,  8,  16,  6,  17]. 


1.2  Modeling  of  ICP  facility 

Inductively  Coupled  Plasma  (ICP)  wind  tunnels  are  usually  hard  to  simulate  due  to  the  com¬ 
plexity  of  the  problem.  The  gas  is  injected  into  the  wind  tunnel  and  then  heated  by  means  of 
an  external  inductor,  then  the  flow  reaches  very  high  temperatures  in  the  torch  (~  lOOOOiF). 
That  implies  that  there  are  several  subsystems  interacting  in  this  kind  of  problems:  on  one 
hand,  the  flow  field,  governed  by  the  Navier-Stokes  equations  at  Local  Thermodynamic  Equi¬ 
librium  (LTE)  (1);  on  the  other  hand,  the  electromagnetic  field  created  by  the  inductor.  Both 
subsystems  are  coupled  and  have  to  be  solved  simultaneously. 
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1.2.1  Electromagnetic  field  equations 


The  electromagnetic  field  equations  are  quickly  presented  in  this  section.  The  equation  of  the 
electric  field  induced  by  the  torch  reads: 

d2E  Id  f  dE\  E  „  ^ 

~dx2  +  rdr  \  ~dr  )  ~  r2  ~  aE  =  2^  5  (r  “  r '0  (20) 

'  '  2=1 

Where  the  electric  field  E  is  a  complex  variable.  We  can  take  into  account  that  the  total  field 
is  equal  to  the  sum  of  the  electric  field  induced  by  the  coils  Ev  and  the  electric  field  produced 
by  the  currents  running  into  the  plasma  Ep: 


E  —  Ev  +  Ep 

Thus,  Equation  (20)  can  be  separated  into  two  different  equations: 

(21) 


(22) 

Equation  (21)  can  be  solved  analytically,  whereas  the  solution  of  equation  (22)  is  more  difficult 
to  be  obtained,  and  one  needs  an  electromagnetic  field  model  to  compute  Ep. 

1.2.2  Coupled  system 

The  Navier-Stokes  equations  and  the  electromagnetic  field  equations  are  coupled  and  solved 
together.  The  torch  creates  an  electric  field  that  interacts  with  the  ionized  plasma,  affecting 
its  movement  by  means  of  the  Lorenz  force.  Besides,  the  electromagnetic  field  produces  also 
changes  in  the  flow,  increasing  its  temperature  because  of  the  Joule  effect.  On  the  other  hand, 
the  temperature  changes  influence  the  electric  conductivity  of  the  air,  meaning  that,  in  equation 
(20),  the  induced  electric  field  will  be  also  affected  by  the  temperature  changes.  To  sum  up, 
these  two  systems  of  equations  are  coupled  in  both  directions:  the  solution  of  the  electric  field 
system  affects  the  Navier-Stokes  system  and  vice  versa. 

A  detailed  explanation  about  the  computation  and  the  implementation  of  the  Lorentz  force 
and  the  Joule  effect  can  be  found  in  [18]  and  is  detailed  bellow. 


1  d 


OEr 


dr  \  dr 
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d*Ev 


Er, 


1  d  (  dE% 


dr  \  dr 


dx 2 
d2En 


-iu)HoIc  E^ 


2=1 


dx 2 


Ep 

=  iunocr  (. Ev  +  Ep 


Lorentz  force.  The  motion  of  a  particle  is  governed  by  the  Lorentz  force,  which  acts  on  a 
point  charge  due  to  the  electromagnetic  field.  It  is  given  by  the  following  equation: 

°—=FL=q{E  +  vxB)  (23) 

where  Fl  is  the  force  (in  newtons),  E  and  B  are  the  electric  (in  volts  per  meter)  and  magnetic 
field  (in  Tesla’s) ,  q  is  the  electric  charge  of  the  particle  (in  coulombs)  and  v  is  the  instantaneous 
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velocity  of  the  particle  (in  meters  per  second),  mv  the  momentum  of  the  particle.  Combining 
the  following  expression: 

E  =E  +  uxB  (24) 

of  the  electric  field  in  a  frame  of  reference  moving  with  the  mean  velocity  of  the  flow,  with  the 
definition  of  Vi  ,  we  can  write  the  Lorenz  force  applied  to  a  particle  i  as  follows: 

Fl,i  =  qi{E  +  v  x  B)  =  qi(E+(u  +  Vi)  x  B)  ==  ®(£  +  «x  B  +  Vi  x  B)  =  qi(E  +  Vi  x  B)  (25) 

Joule  Heating.  The  heat  source  term  due  to  Joule  effect,  also  known  as  ohmic  heating, 
applies  to  any  conductor  when  electric  current  is  passing  through  it.  It  is  defined  by  the  following 
equation: 

Pj  =  ^  UiQi  <  Vi  >  xB  =  Jc-  Ei ,  (26) 

i 

where  i  is  a  charged  particle,  JT  riiqi  <  F{  •  Vi  >  is  the  power  developed  by  a  force  on  a  particle 
Fi  ,  and  Jc  is  the  conduction  current  density.  This  equation  further  simplifies  to 

Pj  =  oEt  •  Ei  (27) 

The  neutral  plasma  hypothesis  JT  mqi  =  0,  the  Reynolds  magnetic  number  and  the  ambipolar 
diffusion  hypothesis  has  been  used  to  derive  these  expressions.  The  root  mean  square  of  Pj  will 
be  used  as  source  term  in  the  N-S  system,  and  it  can  be  easily  derived 

<  Pj  >=<  EI-EI>=a<Ee-Ee  =  a<  E2ec(elujt)22  >=  -^E2^.  (28) 

The  implemented  equation  uses  the  real  and  imaginary  part,  as  follows: 

<  Pj  >=  2  a(pRe  +  P'6, elm)  (29) 


2  Numerical  method 

2.1  Space  discretization 

All  the  results  to  be  presented  have  been  obtained  by  means  of  a  parallel  FV  solver  for  unstruc¬ 
tured  grids  implemented  within  COOLFluiD.  The  FV  discretization  is  applied  to  the  system 
of  governing  equations  written  in  integral  conservation  form: 

f  U  dn  +  <f  Fc  •  n  ddn  =  <f  Fd  ■  n  ddtt  =  [  S  dti  (30) 

ot  Jn  JddQ  JddQ  Jn 

where  Fc  and  represent  respectively  the  convective  and  diffusive  fluxes,  while  S  contains  the 
reaction  terms.  We  apply  a  conventional  cell  centered  approximation,  which  assumes  solution 
vectors  located  at  the  centroid  of  each  computational  cell.  Inverse-distance  weighted  least  square 
reconstruction  [19]  is  utilized  to  yield  second  order  accuracy.  Oscillation  free  solutions  are 
obtained  with  the  multidimensional  limiter  of  Venkatakhrisnan  [20] .  A  modified  version  of  the 
AUSM+  scheme  [21,  4]  has  been  used  to  discretize  the  convective  fluxes  in  all  our  computations, 
as  detailed  hereafter. 
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2.1.1  AUSM+  scheme 


The  AUSM+  scheme  [21]  offers  a  good  compromise  between  robustness  and  accuracy,  while 
being  not  prone  to  the  carbuncle  phenomenon.  The  AUSM+  scheme  relies  on  the  splitting 
between  a  convective  and  a  pressure  component: 

F1/2(Ul,  Ufl,  n)  =  F<c>  +  Fto  =  m1/2*L/R  +  p1/2  (31) 

Herein,  in  the  case  of  a  gas  mixture  in  thermal  and  chemical  nonequilibrium  with  Ns  chemical 
components  and  Nm  molecules,  which  is  the  case  of  interest  here,  the  scalar  mass  flux  m,  the 
interface  Mach  number  Mi/ 2,  the  vector  quantity  \l/  and  the  pressure  flux  can  be  expressed 
respectively  as: 


rhi/2  —  Mi/2CLi/2 


PL 

PR 


if  Mi  1 2  >  1, 
otherwise 


(32) 


Mi/2  =  M+(Ml )  +  M~(Mr)  (33) 

*  =  [ya>  v,  H,  ymEvm\  (34) 

P1/2  =  V+{ML)pLn  +  V~(MR)pRn  (35) 

where  the  actual  definition  of  the  split  Mach  number  AV  and  the  split  pressure  functions  V± 
can  be  found  in  [21].  The  interface  speed  of  sound  aq/ 2  appearing  in  Eq.  32  is  defined 


ai/2  =  min  (aL,aR), 


a  — 


ma x(a*,  |  qn  |) 


(36) 


where  for  the  critical  speed  of  sound,  originally  defined  as  a*  =  y  2^+1^  H  in  [21],  we  adopt  a 
more  suitable  definition  [4]  in  our  computations  dealing  with  flows  in  thermo-chemical  nonequi¬ 
librium: 


* 

a  — 


i  27(7  - 1)  H 

27 + 7(7  - 1) 


(37) 


where  the  frozen  specific  heat  ratio  7 
7  =  1  +p/pe  are  considered. 


Ys^ysdhs/dT 

Usdes/dT 


and  the  equivalent  specific  heat  ratio 


2.1.2  Discretization  of  diffusive  Fluxes 

The  discretization  of  the  diffusive  term  in  Eq.  30  leads  to: 

A Tf 

i  Fd-ndE  =  y;G fVf  (38) 

Jt,  /=1 

where  =  F^  -  and  the  diffusive  fluxes  Fd  typically  depend  on  the  variables  P  and  their 
gradients,  i.e.  Fd  =  Fd(P,  VP). 
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Gradient  calculation:  diamond  control  volume.  The  application  of  Green-Gauss’  theo¬ 
rem  within  a  chosen  control  volume  can  be  used  to  determine  the  above  mentioned  gradients: 


VP  = 


hL 


vp  dnv  = 


P  ndEv 


E” 


(39) 


A  popular  choice  for  iiu  on  unstructured  meshes  is  a  diamond-shaped  volume  [22,  23]  like  the 
one  in  Fig.  12  which  is  built  around  the  considered  face  and  which  includes  all  the  face  nodes, 
the  left  and  right  cell  centers  as  vertexes. 


Figure  12:  Diamond  control  volume  for  the  calculation  of  the  gradients  for  the  diffusive  fluxes. 


The  discretized  version  of  Eq.  39  becomes: 


VP 


1 


Nf 


f= 1 


with  P  f  — 


nI 


EE 


3= 1 


(40) 


—  -P 

where  P/  is  a  face-averaged  value  of  P  calculated  from  the  values  Pj  in  the  vertices  of  the 
diamond  volume,  whose  number  of  faces  is  Nf.  Moreover,  Nh  is  the  number  of  vertexes  in  each 
face  /  of  the  control  volume,  i.e.  Nl  is  equal  to  the  space  dimension  of  the  problem  (2  or  3). 


Gradient  calculation:  Kim’s  method.  Alternatively,  the  method  extensively  described 
in  [24]  can  be  used  to  compute  the  gradients.  In  this  case,  the  cell-based  solution  gradients, 
computed  via  a  least  square  reconstruction,  are  corrected  and  projected  onto  the  faces.  This 
method  is  known  to  provide  superior  accuracy  and  robustness  on  high-aspect  ratio  cells. 


2.1.3  Discretization  of  source  Terms 

The  discretization  of  the  source  term  appearing  on  the  left-hand-side  of  30  is  based  on  the  cell 
centered  value  in  a  given  cell  i: 


s(p)  dn  «  s(p i)Qi  =  s a 


(41) 


When  the  source  term  includes  derivatives  of  some  dependent  variable  p  (e.g.  the  stress  term 
tqq  appearing  in  the  case  of  axisymmetric  Navier-Stokes),  these  are  calculated  by  applying  the 
Green-Gauss  theorem,  similarly  to  what  explained  in  Sec.  2.1.2: 
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(42) 


=  ~~  f  Vp  dQi  =  /  p  n  d^i 

*Li  JQi 

In  this  case,  however,  the  chosen  control  volume  coincides  with  the  volume  of  the  current  cell, 
Qi,  while  a  diamond-shaped  one  was  used  for  the  computation  of  the  diffusive  fluxes.  The 
discretized  version  of  Eq.  42  is  identical  to  Eq.  40. 

2.1.4  Boundary  conditions  for  nonequilibrium  solver 

Radiative  equilibrium  wall.  The  radiative  equilibrium  wall  boundary  condition  consists  in 
imposing  that  the  heat  released  from  the  gas  into  the  wall  by  conduction  ( qc°nd )  and  convection 
(< qc°nv )  is  exactly  balanced  by  the  heat  lost  by  radiation  (qr)  from  the  wall  itself.  This  translates 
into  the  non  linear  equation: 


Q(TW)  =  q™nd  +  q™nv  +  qr  =  0  (43) 

to  be  solved  iteratively  at  the  wall.  After  having  substituted  the  actual  expressions  for  all  terms, 
this  expression  becomes: 


Q{TW)  = 


\dT  ^  dTZ' 
dn  ^  m  dn 

\  m  , 


Ns 

+  ^  ha  Js  •  nw  -  o {eewT^ 


S 


(44) 


where  a  =  5.67 •  10-8  [W/m2 / KA]  is  the  Stefan-Boltzmann  constant,  tew  and  are  respectively 
the  wall  emissivity  and  absorptivity  and  I ^  is  the  distant  body  temperature.  The  solution  Tw 
of  the  non  linear  equation  Q(TW )  =  0  is  obtained  by  applying  a  Newton  procedure,  where  in 
order  to  enhance  numerical  robustness  the  maximum  variation  of  wall  temperature  for  each 
boundary  face  between  two  subsequent  time  steps  is  limited  to  a  user-defined  value  (typically 
100-200  K). 


Super  catalytic  wall  with  LTE.  This  condition  imposes  the  chemical  composition  of  the 
gas  mixture  at  the  wall  to  be  equal  to  the  equilibrium  values.  Under  the  local  thermodynamic 
and  chemical  equilibrium  assumption,  the  species  molar  composition  at  the  wall  Xw  can  be 
computed  from  the  local  temperature  Tw,  pressure  pw  and  elemental  molar  composition  X!^ 
by  solving  iteratively  an  algebraic  system  of  equations  [4,  25,  26].  In  our  implementation,  pw  is 
extrapolated  from  the  interior  cell,  whilst  the  wall  temperature  is  calculated  from  the  radiative 
equilibrium  condition,  as  explained  in  the  previous  section.  While  in  a  full  non-catalytic  case, 
the  heating  due  to  gas  convection  is  null,  i.e.  qg°nv  =  hs^s  •  =  0,  if  we  assume  local 

equilibrium  at  the  wall,  this  term  plays  a  role,  since  the  gradient  of  species  fractions  across  the 
wall  is  not  null  and  contributes  to  increase  the  heat  flux  similarly  to  a  full  catalytic  condition. 


Catalytic  wall  with  finite  rate  chemistry.  A  catalytic  wall  promotes  recombination  of 
atoms  though  exothermic  reactions.  For  this  reason  it  is  very  important  to  account  for  this 
phenomenon.  On  a  partially  catalytic  wall  we  need  to  impose  the  molar  fraction  of  every 
species  i  such  that  the  balance  between  the  diffusive  flux  and  the  production  rate  is  respected. 
To  do  this  we  solve  the  following  equation  by  a  Newton  Method: 

Ji  (xj=i,Ne )  =  Ui  (xj=ltNs )  (45) 
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where  Ji  represents  the  diffusion  flux  of  the  spices  i,  (Ji  the  production  rate,  and  Ns  is  the 
number  of  species  considered.  The  diffusive  fluxes  are  computed  using  Fick’s  law.  In  this  case, 
the  diffusive  flux  Ji  is  proportional  to  the  spatial  derivative  of  the  molar  fraction  of  the 
species  i : 

Ji  =  pDFick ^  (46) 

In  order  to  ensure  the  conservation  of  mass,  the  diffusive  coefficient  DFlck  of  Fick’s  law  is  the 
same  for  all  species  and  it  is  defined  by: 


j-^Fick  _ 


i  -m 


s  = 


Na 


js= o  "'3' 


(47) 


where  yi  is  the  mass  fraction  of  an  arbitrary  species  and  the  Dij  are  the  binary  diffusion 
coefficients.  The  wall  rate  production  is  defined  by: 

F  =  Mjji  r)  -  J2  E  i,  r )  (48) 

r  j  r 


where  (7z)zg[o,ws]  is  the  catalytic  recombination  factor  that  is  identified  by  the  experiments. 
M.\  is  the  impinging  flux  of  species  i  on  the  wall  and  it  is  defined  by: 


M\  =  yip 


TWR 
211  rrij 


(49) 


where  Tw  is  the  temperature  of  the  wall  and  mi  is  the  molar  mass  of  species  i.  Finally,  matrices 
v  and  /I  define  the  recombination  reactions  that  happen  at  the  wall: 


•  z/(i,  r)  =  1  if  the  species  i  is  destroyed  during  the  reaction  r 

•  /i(j,  i,r)  =  1  when  species  j  is  destroyed  during  reaction  r  it  produces  species  i 


2.1.5  Boundary  conditions  for  ICP  solver 

Symmetry  axis.  The  electric  field  has  to  vanish  on  the  symmetry  axis: 

Ep,Re  Ep,Im  =  0.  (50) 


2D  boundary.  We  use  the  following  formula  to  compute  Ep 

cells 


T7I  •  P’0  \  ^ 

E 0  =  "2 *  2. 


2=1 


—jiSiG(ki 

r 


(51) 


where  u  is  the  angular  frequency,  /io  is  the  magnetic  permeability  of  free  space  ;  (r,  z)  are  the 
coordinates  of  the  point  where  the  Electric  Field  is  to  be  computed;  the  summation  extend 
over  the  number  of  cells  i,  (77  ,  Z{  )  are  the  coordinates  of  the  ith  element  of  the  domain  and 
Si  its  area,  and  ji  is  the  current  density  in  it.  And  G  is  defined  by: 


G(k)  = 


(2  -  k2)K(k )  -  2E(k) 


k 


k  = 


4. Rcr 


(■ Rc  +  r )2  +  (z  -  Zc) 2 


Where  K(k )  and  E(k)  are  the  complete  elliptic  integrals  of  the  first  and  second  kind. 
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2.1.6  Implicit  time  integration 


An  implicit  Backward  Euler  method  is  applied  to  Eqs.  (30)  in  order  to  converge  to  steady  state: 


EE1  +  uFV(p)  =  o  =  r(p) 

where  R(P)  is  a  pseudo-steady  residual,  U  =  U(P)  an  explicit  analytical  relation  and  the 
vector  of  natural  variables  P  =  [ps,  v,  T,  Ty],  in  which  the  governing  equations  are  explicitly 
closed  and  in  which  it  is  therefore  convenient  to  store  the  solution  vector.  The  application  of  a 
one  step  Newton  method  yields  the  following  linear  system: 

APn  =  -R(Pn),  (53) 

where  the  jacobian  matrix  ^  is  computed  numerically.  The  GMRES  algorithm  with  an  Addi¬ 
tive  Schwartz  preconditioner  serves  to  solve  in  parallel  the  corresponding  linear  systems  arising 
from  Newton  linearizations.  The  solution  update  is  also  performed  in  primary  variables. 


3  Numerical  results 

3.1  ICP  simulation 

The  first  step  of  the  project  was  to  simulate  the  flow  in  the  Plasmatron  in  order  to  provide  the 
inflow  data  at  the  exit  of  the  torch.  The  chosen  test  conditions  correspond  to  a  static  pressure 
Ps  —  2010  and  a  plasma  power  of  80 kW .  The  boundary  conditions  are  shown  in  Fig.  13. 


subsonic  outlet 

/ 


Figure  13:  Setup  of  the  boundary  conditions  for  the  ICP  simulation. 
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Figure  14:  Temperature  isolines  (top)  and  stream  lines  (bottom)  for  the  ICP  simulation. 
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3.2  Nonequilibrium  simulation 


The  focus  of  this  project  was  to  simulate  a  VKI  Plasmatron  experiment,  including  nonequilib¬ 
rium  and  gas-surface  interactions  effects.  The  test  consists  in  a  plasma  jet  coming  out  from  a 
ICP  torch  and  impinging  on  a  probe  featuring  different  TPS  materials.  The  actual  TPS  sample 
before  the  testing  is  shown  in  Fig.  15  and  includes  three  patches  with  two  different  catalycities 
(low  cat  -  high  cat  -  low  cat).  The  actual  experiment  is  shown  in  Fig.  16:  the  test  model,  a 
flat  plate  probe  with  a  blunted  leading  edge,  is  schematically  depicted  in  Fig.  16(a),  while 
Fig.  16(b)  provides  a  infrared  (IR)  visualization  of  the  temperature  on  the  probe  during  the 
testing  (courtesy  from  [27]).  For  our  numerical  investigation,  the  geometry  has  been  assumed 
to  be  2D,  even  though  traditional  effects  cannot  be  excluded  in  the  real  case.  A  global  and  a 
detailed  view  of  the  2D  computational  mesh  used  for  this  case  are  shown  in  Fig.  17(a)  and 
Fig.  17(b)  respectively.  A  refinement  in  the  streamwise  direction  has  been  applied  near  the 
junctions  between  different  materials  in  order  to  better  capture  the  jumps  in  temperature  and 
heating.  A  first  cell  size  of  10-5  [m]  has  been  imposed  close  to  the  wall. 


Figure  15:  Three-patches  TPS  sample  with  two  different  catalycities  [27]. 


(a)  Flat  plate  probe  configuration  for  TPS  off- 
stagnation  testing. 


(b)  IR  visualization  during  actual  testing  on  the  flat 
plate  probe. 


Figure  16:  Flat  plate  probe  configuration  and  IR  visualization  during  testing  [27]. 
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An  axisymmetric  ICP  simulation  of  the  full  experiment  (torch  +  chamber  +  test  model)  un¬ 
der  a  simplified  LTE  assumption  (using  the  model  detailed  in  Section  1.2)  has  provided  the 
input  jet  profile  for  the  nonequilibrium  2D  calculation.  In  order  to  initialize  the  latter,  the  flow 
properties  of  the  inlet  subsonic  plasma  jet  (in  LTE  conditions)  have  been  extrapolated  to  all 
domain  in  the  x  direction.  Herein,  in  order  to  stabilize  the  computation,  an  artificial  blowing 
velocity  ux  —  50  [m/s]  has  been  imposed  in  the  inlet  profile  outside  the  jet,  while  uy  =  0  [m/s] 
has  been  strongly  set  on  the  whole  inlet  boundary. 

Following  a  conservative  (probably  even  over-conservative)  approach,  the  nonequilibrium  cal¬ 
culation  has  been  run  in  three  steps  using  a  limited  second  order  scheme  from  the  start: 

1.  At  first,  a  fully  adiabatic  solution  has  been  obtained  starting  from  the  previously  described 
initial  conditions. 

2.  Then,  a  solution  with  fixed  temperature  and  radiative  equilibrium  on  different  wall 
patches  has  been  computed. 

3.  Finally,  the  desired  catalytic  boundary  conditions  have  been  applied  to  the  different  wall 
patches,  i.e.  super-catalycity  on  the  isothermal  wall  and  the  chosen  partial-catalycity 
profile  in  combination  with  radiative  equilibrium  on  the  rest  of  the  probe. 

In  the  following,  only  the  final  solution  is  presented. 


x  [m]  x  [m] 


(a)  Global  view  of  the  2D  computational  mesh.  (b)  Zoom  near  the  blunted  leading  edge  of  the  plate. 

Figure  17:  Views  of  the  mesh  (20450  cells)  for  the  blunted  plate. 

A  global  view  on  the  computed  flowfield  is  presented  in  terms  of  temperatures  and  electron 
mass  fraction.  Roto-translational  and  (superimposed)  vibrational  temperature  fields  in  Fig.  18 
appear  to  be  rather  similar,  which  suggests  that  thermal  nonequilibrium  effects  are  playing  a 
minor  role  in  the  flow.  In  particular,  the  flow  appears  to  be  in  equilibrium  close  to  the  inlet 
boundary,  which  verifies  a-posteriori  the  LTE  assumption  for  the  flow  coming  out  from  the 
plasma  torch.  Figure  19  shows  that  only  the  core  of  the  jet,  where  temperatures  reach  up  to 
8000  [K],  is  characterized  by  ionized  flow. 
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Figure  18:  Roto-translational  (contours,  solid  isolines)  and  superimposed  vibrational  tempera¬ 
ture  (dash-dot  isolines). 
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Figure  19:  Electron  mass  fraction  field. 

Figure  20  shows  a  comparison  between  numerical  results  computed  with  the  two  gradient 
calculation  methods  described  in  Section  2.1.2  in  order  to  assess  their  influence  on  the  surface 
temperature  and  heat  flux  distributions.  Kim’s  method  (indicated  as  grad  model  1 )  tends  to 
predict  up  to  3%  higher  temperature  where  radiative  equilibrium  is  imposed  at  the  wall  and 
up  to  7%  higher  heat  flux  along  almost  all  the  body.  In  particular,  Kim’s  method  predicts  1.19 
versus  1.11  [MW/m2]  peak  heating  on  the  blunted  leading  edge  and  0.87  versus  0.82  [MW/m2] 
on  the  junction  between  low-  and  high- cat aly city  patches.  The  results  here  presented  show  a 
qualitative  and  quantitative  behavior  (particularly,  in  terms  of  heat  flux  distribution  and  peak 
values)  consistent  with  the  experimental  measurements  in  [27],  even  though  the  conditions 
for  the  experimental  setup  do  not  fully  match  the  simulated  conditions.  A  more  quantitative 
validation  of  the  numerical  results  against  available  experimental  data  will  be  addressed  in  the 
future,  after  the  inclusion  of  a  more  accurate  diffusion  model  for  characterizing  wall  catalycity. 
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Figure  20:  Surface  distribution  of  temperature  and  heat  flux  along  the  plate  computed  with 
two  different  methods  for  gradient  calculation:  Kim’s  (grad  model  1)  and  diamond  (grad  model 
2). 
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Part  III 


Numerical  simulation  of  plasma  jet  around 
the  VKI  rounded  nose  flat  plate  probe 

1  Introduction 

The  aim  of  this  contribution  has  been  to  compute  the  flowfield  around  the  Von  Karman  Institute 
(VKI)  rounded  nose  flat  plate  probe.  The  probe  is  conceived  to  simulate  a  vehicle  off-stagnation 
region  flowfield  inside  the  VKI  Plasmatron  facility.  The  probe,  shown  in  Fig.  21,  consists  of 
a  copper  plate  with  a  12.5  mm  nose  radius  and  can  host  a  152x30  mm  rectangular  insert  at 
26  mm  distance  from  the  leading  edge  [28].  The  insert  can  be  made  of  different  materials;  in 
Fig.  22  the  low-high-low  catalycity  insert  is  shown.  The  high  catalycity  patch  is  realized  by 
coating  the  sample  with  a  Mullite  section  starting  at  50  mm  from  the  insert  upstream  edge 
and  ending  at  110  mm,  the  low  catalycity  part  is  made  of  SPS  C/SiC  [28]. 

CMC  Sample 


Figure  21:  Flat  plate  probe  configuration  for  off-stagnation  testing. 


Figure  22:  Low- high-low  catalycity  insert. 

The  experiment  for  the  flat  plate  probe  has  been  carried  out  for  a  Plasmatron  pressure  chamber 
of  1500  Pa.  The  estimated  enthalpy  of  the  high  temperature  jet  impinging  on  the  probe  is  of 
19.8  MJ/kg  [28].  The  experimental  values  of  surface  temperature  over  the  catalytic  insert  are 
shown  in  Fig.  23.  The  reported  values  are  for  three  different  insert,  a  low  catalytic  one,  a 
high  catalytic  one  and  the  low-high-low  catalytic  one  of  Fig. 22.  In  Sec.  3  we  will  show  the 
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computations  performed,  for  the  low-high-low  insert,  with  the  finite  volume  code  Cosmic  [29]. 


Figure  23:  Experimentally  measured  temperature  over  catalytic  insert. 


2  Numerical  model 

The  computations  of  the  flow-field  around  the  flat  plate  probe  with  the  low-high-low  catalytic 
insert  have  been  performed  with  the  code  Cosmic  [29] .  Cosmic  is  a  multiblock  structured  finite 
volume  code  able  to  deal  with  arbitrary  mixtures  of  chemically  reacting  gases.  The  convective 
fluxes  are  discretized  with  the  artificially  upstream  flux  vector  splitting  (AUFS)  scheme  [30]. 
The  scheme  has  been  preconditioned  by  the  author  [31]  and  is  able  to  compute  flows  over  a  wide 
Mach  number  range,  from  practically  incompressible  flow  to  hypersonic  one.  Diffusive  fluxes 
(stress  tensor,  heat  flux  and  species  diffusion)  are  discretized  with  a  central  scheme  [29].  The 
mixture  thermodynamic  and  transport  properties  are  computed  with  the  PEGASE  library  [32] . 
The  species  diffusion  fluxes  are  computed  with  the  Stefan-Maxwell  equations  in  order  to  ensure 
a  correct  computation  of  the  heat  flux  load  [33]. 

2.1  Boundary  conditions 

The  flat  plate  probe  wall  is  made,  in  the  case  under  study,  by  three  different  materials:  water 
cooled  copper,  the  low  catalycity  material  and  the  high  catalycity  material.  In  the  latter  two 
cases  we  can  assume  that  the  wall  reaches  radiative  equilibrium. 

Over  the  copper  section  the  boundary  condition  for  the  energy  equation  is  very  simple:  the 
wall  temperature  is  assumed  to  be  constant,  because  of  the  water  cooling.  Over  the  insert  we 
assume  radiative  equilibrium  and  impose  that  the  heat  flux  from  the  gas  to  the  wall  is  balanced 
by  the  radiative  heat  flux  lost  by  the  wall. 

Ns 

-AVT  •  nw  +  ^  hiJi  ■  nw  =  ae(T. £  -  T^) 

2=1 
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(54) 


Here  A  is  the  gas  thermal  conductivity,  Ns  the  number  of  species,  vecJi  the  i  —  th  species 
diffusion  flux,  a  the  Stefan-Boltzmann  constant,  e  the  wall  emissivity,  nw  the  wall  normal; 
here  we  have  implicitly  assumed  that  wall  absorptivity  is  equal  to  wall  emissivity.  Eq.  54  is  a 
non-linear  equation  that  can  be  solved  to  determine  the  wall  temperature  Tw. 

The  species  diffusion  fluxes  at  the  wall  are  themselves  computed  by  the  catalycity  boundary 
condition:  the  net  amount  of  species  i  produced  or  consumed  by  surface  reactions  (wiiCat)  has 
to  be  balanced  by  the  diffusion  flux  of  species  i  itself.  Therefore  the  boundary  condition  is: 

Ji,w  *  =  ^ i,cat  (55) 


The  wall  reaction  rate  w^cat  is  computed  as: 

Tb'P  nr 

^i^cat  —  TTiiNl^  ^  ^  V ki^lk  ~  EE  (Hij'vi  rrijM j  (56) 

k= l  l=i  j= 1 


Here  nr  is  the  total  number  of  wall  reactions  and  the  recombination  probability  of  every  wall 
reaction.  The  first  term  in  the  latter  reaction  rate  expression  describes  the  depletion  of  species 
i  due  to  wall  chemical  reactions;  the  second  term  describes  the  creation  of  species  i  by  wall 
reactions.  The  matrix  has  entry  (£;,  i )  equal  to  one  if  the  ith  species  is  destroyed  in  the  kth 
reaction,  otherwise  zero;  the  /q  matrix  has  entry  (z,  j)  equal  to  one  if  the  jth  species  produces 
the  ith  one  in  the  Ith  reaction,  otherwise  zero.  Aij  be  the  number  flux  of  species  z  impinging 
the  surface  and  its  expression,  computed  from  the  kinetic  theory  of  gases,  reads: 


M\  =  Hi 


(57) 


Our  implementation  of  the  catalycity  boundary  condition  of  Eq.  55  is  based  on  the  observation 
that,  at  the  wall,  the  diffusion  flux  is  already  known  and  it  is  equal  to  the  heterogeneous 
recombination  rate  Wi,cat-  The  Stefan-Maxwell  equations  are  multiplied  by  the  normal  to  the 
wall  and  Wi:Cat  is  stuck  in  the  equations,  in  this  way  they  become  a  set  of  equations  for  the 
determination  of  the  species  chemical  composition  at  the  wall  instead  of  species  diffusion  flux. 
This  implementation  of  the  boundary  conditions  has  the  advantage  of  being  fully  consistent 
with  the  model  used  to  computed  species  diffusion  fluxes  in  the  bulk  of  the  flow. 


3  Results 

We  performed  a  series  of  computations  of  the  flow  field  around  the  flat  plate  probe:  our  aim 
is  to  compare  the  computed  wall  temperature  over  the  low-high-low  catalytic  insert  with  the 
experimental  values  [28].  The  computational  grid  is  shown  in  Fig.  24;  as  it  can  be  noticed 
the  grid  covers  only  a  part  of  the  Plasmatron  test  chamber.  VKI  researchers  provided  us  the 
relevant  input  conditions  needed  by  the  Cosmic  code  [34]:  essentially  the  values  of  temperature 
and  velocity  on  the  left  boundary  of  the  computational  domain,  the  chemical  composition  is 
assumed  to  be  in  chemical  equilibrium  at  the  inlet  of  the  domain.  This  assumption  is  also  made 
in  the  rebuilding  procedure  for  the  measurement  of  catalytic  properties  in  the  Plasmatron 
facility  [35,  28].  The  input  values  of  temperature  and  horizontal  component  of  velocity  are 
shown  in  Fig.  25. 

Due  to  the  relatively  low  maximum  temperature  inside  the  plasma  jet  (Tmax  ps  5800  K)  we 
used  five  species  air  (A^,  O2,  NO ,  TV,  O )  for  the  computations:  the  ionization  level  is  not  high 
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Figure  24:  Computational  grid. 
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Figure  25:  Input  boundary  conditions  at  130  mm  from  probe  leading  edge. 


enough  to  significantly  affect  the  heat  flux  [29].  The  reaction  rate  data  of  Park  [36]  were  used 
for  the  chemical  reactions  inside  the  bulk  of  the  flow. 

At  the  wall  we  considered  two  possible  catalytic  reactions:  N  +  N  — ►  N2  and  O  +  O  — >  O2;  we 
assumed  that  the  recombination  probability  7  is  the  same  for  both  reactions.  The  values  of  the 
recombination  rate  probabilities  for  the  three  materials  (copper,  low-cat,  high-cat)  are  shown 
in  Table  3.  We  selected  two  different  values  of  7  for  the  low  catalytic  insert  to  take  into  account 
the  experimental  uncertainties  [28]. 

The  wall  temperature  is  imposed  to  be  equal  to  350  K  on  the  copper  cooled  section  of  the 
probe,  but  it  is  computed  from  the  radiative  equilibrium  condition  on  the  catalytic  insert;  the 
wall  emissivity  is  assumed  to  be  equal  to  0.85  for  both  low  catalytic  and  high  catalytic  sections. 
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Figure  26:  Wall  temperature  over  the  catalytic  insert. 


In  Fig.  26  the  computed  wall  temperature  over  the  catalytic  insert  is  shown,  along  with  the 
experimental  results.  First  of  all  we  notice  a  non  negligible  difference  in  the  temperature  of 
the  low  catalytic  insert  between  case  cl  and  case  c2,  the  temperature  of  the  high  catalycity 
region  being  only  marginally  affected  by  the  catalytic  values  of  the  previous  section.  We  remark 
a  decent  qualitative  agreement  with  the  experiments,  although  the  computed  temperature  of 
the  low  catalytic  region  is  appreciable  lower  than  the  experimental  one.  On  the  opposite  the 
computed  and  experimental  temperatures  of  the  high  catalytic  region  show  a  better  agreement. 
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